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ABSTRACT 

We provide additional information on our recent study of the electromagnetic emission produced during the 
inspiral and merger of supermassive black holes when these are immersed in a force-free plasma threaded by 
a uniform magnetic field. As anticipated in a recent letter, our results show that although a dual-jet structure is 
present, the associated luminosity is ~ 100 times smaller than the total one, which is predominantly quadrupo- 
lar. We here discuss the details of our implementation of the equations in which the force-free condition is not 
implemented at a discrete level, but rather obtained via a damping scheme which drives the solution to satisfy 
the correct condition. We show that this is important for a correct and accurate description of the current sheets 
that can develop in the course of the simulation. We also study in greater detail the three-dimensional charge 
distribution produced as a consequence of the inspiral and show that during the inspiral it possesses a complex 
but ordered structure which traces the motion of the two black holes. Finally, we provide quantitative estimates 
of the scaling of the electromagnetic emission with frequency, with the diffused part having a dependence that 
is the same as the gravitational-wave one and that scales as L n °?~ coli ~ Q 10 / 3-8 / 3 wn ile the collimated one 

& EM ' 

scales as « fi 5 / 3 " 6 / 3 , thus with a steeper dependence than previously estimated. We discuss the impact 
of these results on the potential detectability of dual jets from supermassive black holes and the steps necessary 
for more accurate estimates. 



1. INTRODUCTION 

The gravitational interaction among galaxies, most of 
which are supposed to host a supermassive black hole (BH) , 
with M > 1O 6 M dShankar et al |2004| jLou & Jiang|2008b, 
is a well-established observational fact (Gopal-Krishna et al. 
20031 |Ellison et al ||20TT] |Mohamed & Reshetn ikov||2011[ 
Lambas et al. 2 012[ ). Moreover, in a few documented astro- 



physical cases, strong indications exist to believe that a bi 
nary merger among supermassive BHs has occurred or is on 
going (Rodrigu ez et al.|20 06 ; Komo ssa et al.|200"3||Dotti et al. 
2009). 

A strong motivation for studying supermassive binary black 
holes (SMBBHs) comes from the fact that their gravitational 
signal will be detected by the planned Laser Interferometri c 
Space Antenna (eLI SA/NGO; |Amaro-Seoane et al.| ( f2012| ); 
|Binetruy et al.| ( [2012[ )). When combined to the usual electro- 
magnetic (EM) emission, the detection of gravitational waves 
(GW) from these systems will provide a new tool for test 



ing a number of fundamental astrophysical issues (Cornish 
& Porterl2007jjHaiman et al.|2009| |Phinney|[2009| . For this 
reason, SMBBHs are currently attracting a widespread in- 
terest, both from an observational and a theoretical point of 
view (|Rezzolla|2 009 ; Reisswig et al |2009[[Kesden et al 1 20 1 0J 



Kocsis et al.||201 1| |Tanaka et al.||2012[ |Sesana et al.||2012[ 



Barausse 2 0l2| l. According to the simplest picture that has 



gradually emerged through a series of semi-analytical stud- 
ies and numerical simulations (Milosavljec & Phinney 2005, 



MacFadyen & Milosavljevic 2008; Roedig et al. 2011; Bode 
et al.|2012| ), the accretion disk formed around the two merg- 
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ing BHs, commonly referred to as the "circumbinary" accre- 
tion disk, can follow the dynamical evolution of the system up 
until the dynamical timescale for the emission of GWs, which 
scales like ~ D 4 , where D is the separation of the binary, be- 
comes shorter than the viscous timescale, which instead scales 
like ~ D 2 . When this happens, the circumbinary accretion 
disk is essentially decoupled from the binary, which rapidly 
enters the final stages of the inspiral. Under these conditions, 
neglecting the inertia of the accreting fluid can be regarded 
as a very good approximation. In contrast, magnetic fields 
generated by the circumbinary accretion disk could play an 
important role and the dynamics of the plasma in the inner re- 
gion can then be described within the force-free (FF) approxi- 
mation. These physical conditions are indeed similar to those 
considered in the seminal investiga tions of BH electrod ynam- 
ics of Blandford and Znajek (Blandford & Znajek 1977 1, who 
addressed the question of whether the rotational energy of an 
isolated BH can be extracted efficiently by a magnetic field. 
After the fir st two-dimensional investigations of Komissarov 
and Barkov (Komissarov 2004; Komissarov & Barkov 2009), 
the numerical study of BH magnetospheres has now entered a 
mature phase in the context of SMBBHs evolution. 

In an extensive analysis, but still in t he absence of cur- 
rents and charges, i.e. , in electro vacuum, |Mosta et al.] | |2010) 
showed that, even though the EM radiation in the lowest I = 2 
and m = 2 multipole reflects the gravitational one, the energy 
emitted in EM waves is ~ 13 orders of magnitude smaller 
than that emitted in GWs for a reference binary with mass 
M = 10 8 M Q and a magnetic field B = 10 4 G, thus casting 
serious doubts about a direct detection of the two different sig- 
nals. However, a series of more recent numerical simulations 
in which currents and charges are taken into account, have 
suggested the intriguing possibility that a mechanism similar 
to the original one proposed by Blandford and Znajek m ay be 
activated in the case of binaries (|Palenzuela et al.|2009 2010| 
I Palenzuela et al.|2010b]al|Moesta et al.|2012| >; note that Palen-| 
zuela et aT| (2010b a); Mo esta et al.| | |2012) also make use of a 
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FF approximation. In particular, the Blandford-Znajek mech- 
anism is likely to be valid under rather general conditions, 
namely even if stationarity and axisymmetry are relaxed and 
even if a non-spinning BH is simply boosted through a uni- 
form magnetic field. Moreover, for such uniform magnetic 
field, the emitted EM flux shows a high degree of collima- 
tion, making the EM counterpart more easily detectable. A 



less op timistic view has emerged recently in Moesta et al. 
( |2012[ ) (hereafter Paper I), where we have shown, through 
independent calculations in which the EM emission was ex- 
tracted at much larger radii, that the dual-jet structure is in- 
deed present but energetically subdominant with respect to 
the non-collimated and predominantly quadrupolar emission. 
In particular, even if the total luminosity at m erger is ~ 100 
times larger than in Palenzuela et al. (2010b), the energy flux 
is only ^8 — 2 times larger near the jets, thus yielding a col- 
limated luminosity that is ~ 100 times smaller than the total 
one. As a result, Paper I indicated that the detection of the 
dual jets at the merger is difficult if not unlikely. 

Here we provide additional information on the results pre- 
sented in Paper I and discuss the details of our implementa- 
tion of the equations in which the FF condition is obtained 
via a damping scheme which drives the solution to satisfy the 
correct condition. We show that this is important for a correct 
and accurate description of the current sheets that can develop 
in the course of the simulation. We also study in greater de- 
tail the three-dimensional charge distribution produced as a 
consequence of the inspiral and show that during the inspiral 
it has a complex structure tracing the motion of the two BHs. 
Finally, we provide quantitative estimates of the scaling of the 
EM emission with frequency, with the diffused part having a 
dependence that is the same as the GW one and that scales 
as £ non - co11 ~ ^10/3-8/3 wn ii e t h e collimated one scales as 

EM ' 

Lt° n ~ £! 5 / 3 ~ 6 / 3 , thus with a steeper dependence than previ- 



ously estimated by Palenzuela et al.| 

The plan of the paper is the following. In Section 
describe the system of equations considered in our analysis, 
with particular emphasis on the treatment of the FF condition, 
while in Section[3]we discuss the different routes to the calcu- 
lation of the EM radiated quantities. In Sectionfflwe present 
the astrophysical setup of a BH binary merger, while SectionB] 
compares different approaches for the enforcement of the FF 
condition. Section [6] is devoted to the presentation of the re- 
sults, and, in particular, to the computation of the luminosity. 
Finally, Section|7]contains the conclusion of our work and the 
prospects for the detection of an EM counterpart to SMBBHs. 

In the rest of the paper, we set c = G = 1, adopt the stan- 
dard convention for the summation over repeated indices with 
Greek indices running from to 3, Latin indices from 1 to 3, 
and make use of the Lorentz-Heaviside notation for the EM 
quantities, in which all \J Air factors disappear. 

2. EVOLUTION EQUATIONS 

We solve the combined system defined by the Einstein and 
Maxwell equations and model either an isolated rotating BH 
or a BH binary inspiralling in quasi-circular orbits. In both 
cases we assume that there is an external FF magnetic field. 
More specifically, we solve the Einstein equations 
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(1) 



where g^ u , and are the Ricci, the metric, and the 
stress-energy tensors, respectively. In addition, we solve the 



following extended set o f Maxwell equations (Komissarov 
[20071 |Palenzuela et al.|2009) : 



(2) 
(3) 

where F„„ is the Faraday tensor, *F fl „ is its dual, 1^ is the 
four-current, and we have introduced a 3+1 slicing of space- 
time, with being the unit (future oriented) timelike vector 
associated with a generic normal observer to the spatial hy- 
persurfaces. 

The set of Maxwell equations |2]i and |3]l is referred to as 
"extended" because it incorporates the so-call ed divergence- 
cleanin g approach, originally presented in |Dedner et al. 
(2002) in flat spacetime, and which amounts to introducing 
two additional scalar fields, \& and that propagate away the 
deviations of the divergences of the electric and of the mag- 
netic fields from the values prescribed by Maxwell equations. 
Such scalar fields are initialized to zero, but are driven into 
evolution as soon as violations of the EM constraints are pro- 
duced. The total stress-energy tensor is composed of a term 
corresponding to the EM field: 



Tf = F^xF 



(F XK F Xh 



(4) 



and of a term due to matter, . However, because the EM 
field is assumed to be FF, Tjf T^f, and the total stress- 
energy tensor is then assumed to be given entirely by Equa- 
tion namely w T^ v , In the rest of our discussion we 
will use the expression "electrovacuum" to denote the case 
when currents and charges of the Maxwell equations are zero. 



Such a scenario was extensively studied in |Mosta et al. (2010) 
and it will be used here as an important reference. In what fol- 
lows we discuss in more detail our strategy for the solution of 
the Einstein equations and of the Maxwell system in an FF 
regime. 

2.1. The Einstein Equations 

For the solution of the Einstein equations we make use 
of a three-dimensional finite-differencing code that adopts 
a conformal-tra celess "3 + 1" BS SNOK formulation of the 
equations (see Pollney et al. (2007) for the full expressions in 
vacuum and Baiotti et al.! (2008 ) for the case of a spacetime 
with matter). The code is base d on the Cactus Computa- 
tional Toolkit (Allen et al. 2000) and employs adaptive mesh- 
refine ment techniques via the Carpet-driver (Schnetter et al. 
2004). For compactness we will not report here the details 
regarding the adopted formulation of the Einstein equations 
a nd the gauge conditions use d, which can however be found 



inl Pollney et^p007l|2OTTV 

We also note that recent developments, such as the use of 
eighth-order finite-difference operators or the adoption of a 
multiblock structure to ext end the size of the wave zon e, have 
been recently presented in Pollney et al. (2009 2011 1. Here, 
however, in order to limit the computational costs and be- 
cause a very high accuracy in the waveforms is not needed, the 
multiblock structure was not used and we have used a fourth- 
order finite-difference operator with a third-order Implicit- 
Explicit Ru nge-Kutta (RKIMEX) integration in time (see 
Section|Z3l). 



2.2. The Maxwell Equations 

The Maxwell equations (|2]i and (3) take a more familiar 
form when expressed in terms of the standard electric and 
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magnetic fields as defined by the following decomposition of 
the Faraday tensor in a 3+1 foliation: 



= n»E v 



i v E» 



B a np , 



n . 



(5) 
(6) 



where the vectors E^ and B^ are purely spatial (i.e. , E^n^ = 
B^n^ = 0) and correspond to the electric and magnetic fields 
measured by the normal (Eulerian) observers. The two ex- 
tra scalar fields $ and <£> introduced in the extended set of 
Maxwell equations lead to two evolution equations for the 
EM constraints, which, we recall, are given by the divergence 
equations 



^^E ^ = q , 
V,B l = , 



(7) 
(8) 



where the electric current has been decomposed in the electric 
charge density q = —n^I^ and the spatial current J; = 
More specifically, these evolution equations describe damped 
wave equations and have the effect of controlling dynamically 
the possible growth of the violations of the constraints and of 
propagating them away from the problematic regions of the 
computational domain where they are produced. 
In terms of and B^, the 3 + 1 formulation of Equations 



|2} and ^ becomes ( |Palenzuela et al.|2 010a ) 



V t B l + e ijk Vj (aE k )+a <S> = aKB i 

V t ^ + a W i E i = aq - , 

P t $ + aV i B i = 

V t q + V l {aJ l ) = aKq, 



V t E l - e^Vj (aB k ) + a 7 iJ Vj * = a K E l - a J 1 , (9) 

(10) 

(11) 
(12) 
(13) 

where T> t = (d t — Cp) and Cp is the Lie derivative along 
the shift vector (3 and K is the trace of the extrinsic curva- 
ture. The charge density q can be computed either through the 
evolution equation (\3\ or by inverting the constraint equa- 
tion |7]). For simplicity, we choose the latter approach, which 
ensures that the constraint ( fT2| ) is automatically satisfied if 
= initially and effectively removes the need for the po- 
tential 

Exploiting now that the covariant derivative in the second 
term of Equations ( fl~0| > and (Hi reduces to a partial derivative, 
i.e. , 

e^ k VjB k = e^ k {djB k + Y 1 ^) = & k djB k , (14) 

and using a standard conformal decomposition of the spatial 
3 -metric 



(15) 



we obtain the final expressions for the extended Maxwell 
equations that we actually evolve 



Vt E l - ^ k [ {dj a)~ lck B c + a ( 4 % k d 3 + d 3 j ck ) B c + a % k d J B c ]=aK E 1 -aJ\ (16) 
V t B i + ^ k e 4 * [ (dj a)% k E c + a ( 4 ~ lck d 3 (j> + dj % k ) E c + a % k dj E c }+ae~ 4 * f> V j $ = aKB i , (17) 
V t <f> + aV l B i = -aK$. (18) 



Clearly, the standard Maxwell equations in a curved back- 
ground are recovered for $ = 0, so that the <£> scalar can then 
be considered as the normal-time integral of the standard di- 
vergence constraint (Rj), which propagates at the speed of light 
and is damped during the evolution. 

As mentioned above, the coupling of the Einstein to the 
Maxwell equations takes place via the inclusion of a nonzero 
stress-energy tensor for the EM fields which is built in terms 
of the Faraday tensor as dictated by Equation Q. More 
specifically, the relevant components of the stress-energy ten- 
sor can be obtained in terms of the electric and magnetic 
fields, that is as 



r = n^T"" = — (E 2 + B 2 ), 
S« = -n^ = ^e ijk EiB k , 



S ij ~ ^0 - 4?r 



-E i E j -B i B j + - 7ij {E 2 + B 2 ) 



(19) 
(20) 

(21) 



where E 2 = E k E k and B 2 = n \o k . 
t can be identified with the energy density of the EM field, 
while the energy flux Sj is the Poynting vector. 
As already discussed in the Introduction, we remark again 



that the EM energies that will be considered here are so small 
when compared with the gravitational binding ones that the 
contributions of the stress-energy tensor to the right-hand- 
side of the Einstein equations ([TJ are effectively negligible 
and thus can be set to zero, reducing the computational costs. 
The fully coupl ed set of the Einstein-Maxw ell equations was 
considered in |Palenzuela et af. (2009 2010) and the compar- 
ison with the results obtained here suggests that for the fields 
below < 10 8 G, the use of the test-field approximation is fully 
justified. 

2.3. Numerical Treatment of the Force-free Conditions 

As commented before, within an FF approximation the 
stress-energy tensor is dominated by the EM part and the con- 
tribution coming from the matter can be considered zero. Fol- 
lowing Palenzu ela et al.| ( (2010a[ ), the conservation of energy 
and momentum, V^'J'^ 0, implies that also the Lorentz 
force is negligible, i.e. , 



(22) 



which can also be written equivalently in terms of quantities 
measured by Eulerian observers as 



E k J k -0, 

qE l + e ijk J j B k = 0. 



(23) 
(24) 
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Computing the scalar and vector product of the equations 
above with the magnetic field B l , we obtain 



E k B k = 

r = q 



B 2 



Jb 



B 1 
B 2 



(25) 
(26) 



The first relation ( |25] > implies that the electric and magnetic 
fields are orthogonal, while expression |26) defines the cur- 
rent, whose component parallel to the magnetic field, namely 
Jb = J l Bi, needs to be defined via a suitable Ohm law. 
From the numerical point of view, specific strategies must be 
adopted in order to enforce the FF constraints expressed by 
Equations ( |25) and (26) . In fact, even though such constraints 
are exactly satisfied at time t = 0, there is no guarantee that 
they will remain so during the evolution of the sys t em. 

The approach introduced by Palenzuela et al. ( 2010a) to 
enforce the constraints ( |2"5) and (|26) consists in a modinca- 
tion of the system at the discrete level, by redefining the elec- 
tric field after each timestep in order to remove any compo- 
nent parallel to the magnetic field. In other words, after each 
timestep the newly computed electric field is "cleaned" by im- 
posing the following transformation (Pa lenzuela et al.|2010a I 



E 1 



E l - (E k B k 



B l 
B 2 



(27) 



In addition, the current is computed from Equation ( |26) af- 
ter settingJe = . An alterna tive approach, introduced 
Komissarov (2011 1 and then in |Lyutikov (2011 1, uses the 



Maxwell equations to compute V t (E B k ), which has to van- 
ish according to Equation ( |25) . Using Equations ( (TO) and ( fTT) 
it is then easy to obtain the following prescription for Jb- 

Jb = \ [B^V^aBu) - Eie iik Vi{aE h )] . (28) 

Without further modifications, however, this approach leads 
to large violations of the FF constraint ( |2"5) in long-term nu- 
merical simulations, as it does not provide a mechanism for 
imposing the constraint at later times. 

As we will show later on, both approaches ( |2"7) and ( |2"8) 
are not fully satisfactory and, as a consequence, we here 
present an alternative method, which takes inspiration from 
the treatment of currents (and related stiff source terms) in 
resistive magnetohydrodynamics. T he idea of in troducing a 
suitabl e Ohm law was proposed in Komissarov ( 2004) and 
then in Palenz uela et al.| ( |2010a) , but it has not been used so 
far in numerical simulations, due to the presence of stiff terms 
which appear as a result. In practice, our continuum approach 
is equivalent to the insertion of suitable driver terms, so that 
the parallel component Jb is computed from an Ohm law of 
the type 



Jb — &BE k Bk, 



(29) 



where as is the anisotropic conductivity along the magnetic- 
field lines. This additional term in the current acts like a 
damping term in the evolution dt(E k B k ), and enforces the 
constraint ( pB") on a timescale l/crg. For <tb sufficiently 
large, one can ensure that the FF constraint ( p5) is always sat- 
isfied. In the simulations presented in this paper, we choose 
<jb > 1/At, where At is the timestep on the finest re- 
finement level. The resulting hyperbolic system with stiff 
terms is solved using a third-order RKIMEX time integration 
method with the technical implementation following the one 



discussed in Pa lenzuela et al.| ( [2009) and with additional de- 
tails presented in the Appendix. 

An additional problem in the numerical treatment of the FF 
approach is represented by the development of current sheets, 
namely of regions where the electric field becomes larger than 
the magnetic field, such that the condition 



B 2 - E 2 > 



(30) 



is violated. If this happens, and in the absence of a proper 
Ohm law responsible for the resistive effects, the Alfven wave 
speed becomes com plex and the system of FF equations is no 
longer hyperbolic (Komissarov 2004). Under realistic con- 
ditions, one expects that in these regions an anomalous and 
isotropic resistivity would restore the dominance of the mag- 
netic field. A solution to this problem was proposed in Komis- 
( 2006), where the velocity of the drift current was mod 



sarov 



ified in order to ensure that it is always smaller than the speed 
of light. This leads to the following prescription for the cur- 
rent: 



J ~ q B 2 + E 2 



Jb 



B l 
£2 



(31) 



which should be compared with Equation ( |26) and has the net 
result of underestimating the value of the current. 

An alternative solution to the numerical treatment of current 
sheets consists in a modification of the system again at the 
discrete level (Palen zuela et al.|2010a) . In practice, after each 
timestep a correction is applied "by hand" to the magnitude of 
the electric field in order to keep it smaller than the magnetic 
field, i.e. , 



E' 



E' 



(i-e) + e^ 



(32) 



with 9 = 1 when B 2 — E 2 < and = otherwise. 

Our strategy, however, differs from both the previous ones 
and follows the same philosophy behind the choice of the 
driver defined by Equation p9) . We therefore introduce a sec- 
ond driver in Ohm law, which will act as a damping term for 
the electric field in those cases when E 2 > B 2 . This addi- 
tional term, combined with the prescription for the parallel 
part of the current ( |29) , leads to the following effective Ohm 
law: 



J 1 



£ l ° k E 3 B k 
B 2 



q- 



a B (E k B k 



- a B (B 2 - E 2 )E^ . 

(33) 

Expression ( (33) shows therefore that in normal conditions, 
i.e. , when B 2 ^ E 2 > 0, the last term introduces a very small 
and negative current along the direction of the electric field. 
However, should a violation of the condition ( (30) take place, 
a positive current is introduced, which reduces the strength of 
the electric field and restores the magnetic dominance. 

In Section[5]we will compare the different prescriptions for 
the enforcement of the FF condition and show that, in con- 
trast to recipes ( |27) and ( (32) , our suggestions ( |29) and ( (33) 
yield both and accurate and a smooth distribution of the EM 
currents. 

3. ANALYSIS OF RADIATED QUANTITIES 

The calculation of the EM and gravitational radiation gen- 
erated during the inspiral, merger and ringdown is an impor- 
tant aspect of this work as it allows us to measure the amount 



5 



correlation between the two forms of radiation. We compute 
the gravitational radiation via the Newman-Penrose curvature 
scalars. In practice, we define an orthonormal basis in the 
three-dimensional space (f, 9, 0), with poles along z. Using 
the normal to the slice as timelike vector t, we construct the 
null orthonormal tetrad {l,n,m, m}: 

l=-^=(t + f), n=-^(t-f), m = -^(0 + i(j>) , 

(34) 

with the bar indicating a complex conjugate. Adopting this 
tetrad, we project the Weyl curvature tensor C a /3-yS to ob- 
tain ^4 = C a j3 1 sn a rh l3 ri 1 m 5 , that measures, ideally at 
null infinity, the outgoing gravitational radiation. For the 
EM emission, on the other hand, we use two equivalent ap- 
proaches to cross-validate our measures. The first one uses 
the Newman-Penrose scalars $o (for the ingoing EM radia- 
tion) and $2 (f or the outg oing EM radiation), defined using 
the same tetrad (Teukols ky|1973| l: 



modified. More specifically we rewrite it as 



$ = F^l v m, 



(35) 



By construction, the Newman-Penrose scalars ^4,^0, $2 are 
dependent on the null tetrad |34|, so that truly unambiguous 
scalars are measured only at very large distances from the 
sources, where inertial observers provide preferred choices. 
Any measure of these quantities in the strong-field region is 
therefore subject to ambiguity and risks to produce mislead- 
ing results. As an example, the EM energy flux does not show 
the expected 1/r 2 scaling when $ 2 and $0 are measured a t 
distances of r ~ 20 M, as used in Palenzuela et al. (20 10b|a| l, 
which is instead reached only for r > 100 M. As we will 
show in Section |6j this fact is responsible for significant dif- 
ferences in the estimates of the non-collimated EM emission. 

The use of a uniform magnetic field within the compu- 
tational domain has a number of drawbacks, most notably, 
nonzero initial values of $ 2 , $o- As a result, great care has 
to be taken when measuring the EM radiation. Fortunately, 
we can exploit the linearity in the Maxwell equations to dis- 
tinguish the genuine emission induced by the p resence of 
the BH (s) from the background one. Following |Teukolsky| 
( |1973| ), we compute the total EM luminosity as a surface in- 
tegral across a 2-sphere at a large distance: 



lim — 

r— foo 27T 



(1*2 



din. 



(36) 



which results straightforwardly from the integration of the 
component of EM stress-energy tensor Q along the time- 
like vector n M and the normal direction to the large 2-sphere 
(namely, the flux of the Poynting vector in Equation ( p"9") 
through the 2-sphere). The term $0 in Equation ( po} has been 
maintained (it disappears at null infinity) to account for the 
possible presence of an ingoing component in the radiation at 
finite distances. In particular, Equation ([36) shows that the 
net flux is obtained by adding (with the appropriate sign) the 
respective contributions of the outgoing and ingoing fluxes. 
More specifically, in terms of the complex scalars $ 2 and 
the outgoing net flux is obtained by subtracting the square of 
their respective moduli. In the specific scenario considered 
here, where a nonzero non-radiative component of the mag- 
netic field extends to large distances, expression (136) must be 



£„„ = lim 



1 

00 2n 



Jr 2 (|$ 2 - $ 2 ,b| 2 - |*o - $o,b| 2 ) dn , 



(37) 

where <£>2,b and ^o.b are the values of the background scalars 
induced by the asymptotically uniform magnetic-field solu- 
tion in the time-dependent spacetime produced by the binary 
BHs. Under assumption of a vanishing net ingoing radia- 
tion, i.e. , $0 ~ *o,b and of stationarity of the background 
fiel d, i.e. , ( &2.b ~ ^o.B; expres sion (|37li can also be rewritten 
as ( |Neilsen et al.|201 i||Ruiz et al.|2QrZ 



lim 



1 

2^ 



(|$2-$o| 2 )^- (38) 



Although Equation (38 1 does not represent, at least in a strict 



physical and mathematical sense, a valid expression for the 
emission of EM radiation in generic scenarios, it can provide 
a useful recipe whenever the assumed approximations made 
above are actually fulfilled. In Section [6] we will assess to 
what degree this is the case for the specific scenario and model 
considered here. 

The choice of the background values of the Newman- 
Penrose scalars <E>2.b and <£>o,b plays a crucial role in measur- 
ing correctly the radiative EM emission, since these quantities 
are themselves timedependent and cannot be distinguished, at 
least a priori, from the purely radiative contributions. This 
introduces an ambiguity in the definition of $ 2 ,b and $o,b^ 
which can however be addressed in at least two different 
ways. The first one consists in assuming that the background 
values are given by the initial values, and further neglecting 
their time dependence, namely setting 



$2,B=$ 2 (i = 0) 



$0,: 



* (t = 0) . 



(39) 



Since all the m = multipoles of the Newman-Penrose 
scalars are not radiative, a second way to resolve the ambi- 
guity is to remove those multipole components from the esti- 
mates of the scalars, namely, of defining 



2.B 



($2)<,m=0 . $0,B = ($okm=0 , 



(40) 



where ($2)e,m=o refer to the m = modes of the multi- 
polar decomposition of <5> 2 (t < 8 is sufficient to capture 
most of the background). Note also that because the rn = 
background is essentially time independent (after the initial 
transient), the choice ( |40] i is effectively equivalent to the as- 
sumption that the background is given by the final values of 
the Newman-Penrose scalars as computed in an electrovac- 
uum evolution of the same binary system. While apparently 
different, expression s (|39| ) and ( |40| > lead to very similar es- 
timates (see Section |6.1| l and, more importantly, they have a 
simple interpretation in terms of the corresponding measures 
that they allow. 

The second approach that we have followed for the compu- 
tation of the emitted luminosity is the evaluation of the flux 
of the Poynting vector across a 2-sphere at large distances 
in terms of the more familiar 3+1 fields E l and B % in Equa- 



tion ( 19 1. Of course, also such evaluation is adequate only far 
from trie binary. The purpose of implementing both versions 
of the luminosity calculation, that are conceptually equivalent 
but differ in the technical details, is precisely to quantify the 
error introduced by evaluating the flux at large but finite dis- 
tances via the Newman-Penrose scalars $2 and $o- Also in 
this case, to account for the background non-radiative contri- 
bution due to our choice of uniform magnetic field (and using 
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again the linearity in the Maxwell equations), we need to re- 
move the background values of the EM fields E^, B 3 B . The 
relevant part of the Poynting vector is then computed as 



(41) 



where, consistently with expression (|39|), we set 



E k (t = 0)=0, 



B k {t 



0) ? . (42) 



As we will show in Sections |6.2| and |67T] we have verified 
that the measures of the EM luminosity obtained using Equa- 
tion |39) or Equation ( |4"0] > reproduces well the corresponding 
ones obtained using the Poynting vector in Equation ( |4T| i. 

4. ASTROPHYSICAL SETUP AND INITIAL DATA 

As mentioned in the Introduction, the astrophysical sce- 
nario we have in mind is represented by the merger of su- 
permassive BH binaries resulting from galaxy mergers. More 
specifically, we consider the astrophysical conditions during 
and after the merger of two supermassive BHs, each of which 
is surrounded by an accretion disk. As the merger between 
the two galaxies takes place and the BHs get closer, a sin- 
gle "circumbinary" accretion disk is expected to form, reach- 
ing a stationary accretion phase. During this phase, the bi- 
nary evolves on the dynamical viscous timescale of the 
circumbinary accretion disk, which is regulated by the abil- 
ity of the disk to transport its angular momentum outward 
(either via shear viscosity or magnetically mediated instabili- 
ties). On a much longer radiation-reaction timescale r GW , the 
system looses both energy and angular momentum through 
the emission of GWs, hence progressively reducing the bi- 
nary separation D. As a consequence, for most of the evo- 
lution the disk slowly follows the binary as its orbit shrinks. 
However, because r GW and Td have a very different scaling 



with D, more specifically r G 



D 4 while r d - D , at 



a certain time the timescale tqw becomes smaller than Td. 
When this happens, the disk becomes disconnected from the 
binary, the mass accretion rate reduces substantially and the 
binary performs its final orbit s in an "interior" region which 
is essentially devoid of gas (|Ar mitage & Natarajan 2002 Liu 
et al. 2003 Milosavljec & Phinney 200"5][ This represents 
the astrophysical scenario in which our simple model is then 
built. 

Although poor in gas, the inner region is coupled to the 
circumbinary disk via a large-scale magnetic field, which we 
assume to be anchored to the disk. The inner edge of the 
disk is at a distance of ~ 10 3 M and is effectively outside 
of our computational domain, while the binary separation is 
only of D ~ 10 M. For simplicity, and because a large-scale 
dipolar field will appear as essentially uniform on the orbital 
lenght scale of the binary during the final stages of the inspi- 
ral, we use an initially uniform magnetic within the computa- 
tional domain. More specifically, the initial magnetic field has 
Cartesian components given simply by B % = (0, 0, Bq) with 
Bq M = 10~ 4 in geometric units or Bq ~ 10 8 G for a binary 
with total mass M = 10 8 M o Q Furthermore, because we 
consider the initial conditions to represent a tenuous plasma 
electrically neutral, the charges, electric currents, and the ini- 
tial electric field are all assumed to be zero, i.e. , E l = = q. 

5 Smaller values of the magnetic field would lead to a less accurate esti- 
mates of the EM fields, but have also been considered. No appreciable differ- 
ences have been measured when using a magnetic field Bq M = 10 — 6 . 



We note that although reasonable, the assumption of a 
large-scale uniform magnetic field has a deep impact on the 
results obtained and more realistic magnetic-field topologies 
will be considered in our future work. As mentioned earlier, 
although astrophysically large, the initial magnetic field con- 
sidered here has an associated EM energy which is several or- 
ders of magnitude smaller than the gravitational-field energy 
and can be are treated as a test field. On the other hand, the 
combination of very low densities and strong magnetic fields 
makes the FF approximation rather appropriate for capturing 
the dynamics of the tenuous plasma. 



4.1. Initial Data and Grid Setup 

We construct consiste nt BH initial d ata via the 
method as described in |Ansorg et al 



(2004). 



puncture" 
We consider 



binaries with equal masses but with two different spin config- 
urations: namely, the s binary, in which both BHs are non- 
spinning, and the s 6 binary, in which both BHs have spins 
aligned with the orbital angular momentum. We use these 
two configurations to best isolate the effects due to the binary 
orbital motion from those related to the spins of the two BHs. 



We note that similar initial data were considered by Koppitz 
[eni^ ( [2007) i; |Pollney et aT| ( [2007] l; [Rezzolla et al.j ( |2008c|b|a| i 
but we have recalculated them here using both a higher reso- 



lution and improved initial orbital parameters. More specifi- 
cally, we use post-N ewtonian (PN) evo lutions following the 
scheme outlined in |Husa et al. [ ( |2008[ ), which provides a 
straightforward prescription for initial-data parameters with 
small initial eccentricity, and which can be interpreted as part 
of the process of matching our numerical calculations to the 
inspiral described by the PN approximations. The free param- 
eters of the puncture initial data are then: (1) the puncture co- 
ordinate locations, (2) the puncture bare mass parameters, (3) 
the linear momenta, and (4) the individual spins. The parame- 
ters of the models adopte d in the numerical si mulations can be 
found in Koppit z et a"L] ( |2007| l; [Pollney et al] ( 12007) ; |Rezzolla 
et al. ( 2008c|b|a[ ). In brief, the initial separation is D = SM 



for all of them, where M is the total initial BH massj^jchosen 
as M = 1, while the individual asymptotic initial BH masses 
are Mj = 1/2. In addition, the EM field is initialized to 

B i = (0,0, B ) with B ~ 10" 4 /M - 10 8 (10 8 M Q /M) G 
and E i = 0. 

The numerical grids consist of nine levels of mesh refine- 
ment, with a fine-grid resolution of Ax/M = 0.025. The 
wave-zone grid, in which our wave extraction is carried out, 
has a resolution of Ax /M = 1.6, and extends from r = 24 M 
to r = 180 M. Finally, the outer (coarsest) grid extends 
up to a distance of ~ 820 M in each coordinate direction. 
Shorter, higher-resolution simulations have also been carried 
out to perform consistency checks. Finally, in addition to 
BHs in a binary system, we have also considered spinning 
and non-spinning isolated BHs as testbeds for our implemen- 
tation of the FF condition. In this case, the numerical grids 
consist of seven levels of mesh refinement, with a fine-grid 
resolution of Ax/M = 0.04 and a coarse-grid resolution of 
Ax/M — 2.56, placing the outer boundary at a distance of 
~ 410 M in each coordinate direction. 

5. ACCURATE FORCE-FREE ENFORCEMENT 



As mentioned in Section 2.3 several different approaches 
are possible to enforce the FF conditions (p5| and d26l) in the 



6 Note that the initial ADM mass of the spacetime is not exactly 1 due to 
the binding energy of the BHs. 
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FIG. 1. — Top row: orthogonality condition (left panel) and current-sheet condition (right panel) for a single spinning BH (dimensionless spin parameter 
a = J/M 2 = 0.7), using different prescriptions for the current: fully discrete approach (light-blue solid line), driveri plus discrete2 (red dotted line), 
driveri plus continuum (dark-blue dashed line), driveri plus driven (black long-dashed line). Bottom row: the same as in the top row, but for the 
equal-mass non-spinning binary BH system sq. 



plasma. The i mportant advantage of the discretized approach 
introduced by Palenzuela et al. ( 2010a| ) is that, at least glob- 
ally, it gives the desired result of an FF solution. In fact, since 
this approach acts "by hand" on the EM fields and converts 
them to values which would yield an FF regime, one is guar- 
anteed that the constraints (25) , ( |2"6) , and ( (30) are satisfied. 
However, a potential disadvantage of such approach is also 
that there is no guarantee that the solution that is forced lo- 
cally with the transformations (|27)-((32) is compatible with 
the solutions in their neighborhoods and thus, that it leads to 
a smooth and accurate representation of the EM fields in the 
presence of current sheets]^] As we will show below, this con- 
cern is indeed well grounded, but it can be resolved effectively 
through the "driver" approach proposed here. 

To compare the different FF prescriptions we have consid- 
ered the simpler setup of a single spinning BH as this allows 
us to concentrate on stationary solutions and hence to isolate 
the potential drawbacks of the different prescriptions, which 
in a binary would otherwise be confused with the actual dy- 
namics of the EM fields. Figure [T] reports the time evolution 

7 Indeed, it is a common experience that any local numerical modification 
of the solution, e.g. , in terms of boundary conditions, is likely to be incom- 
patible with the solution in the bulk. 



of the 2-norms of the scalar product E l Bi, i.e. , ||i? l i?i||2 
(left column) and of the fractional 2-norm of (B 2 — E 2 ), i.e. , 
1 - \\B 2 - E 2 \\ 2 /(\\B 2 - E 2 \\ 2 ) t=0 (right column), moni- 
toring possible deviations from the orthogonality condition of 
Equation ( |25) and from the current-sheet condition of Equa- 
tion ( (30) . The top row of Fig.[T] in particular, refers to a single 
spinning BH, while the bottom row has been obtained in the 
case of the non-spinning BH binary sq. 

The different curves correspond to the various combina- 
tions in the specification of the current and in the treatment 
of the FF constraints. In particular, the labels in the legend of 
Figure [TJrefer to the following choices: 

• discret e^ denotes the first step o f the "discrete" ap- 
proach of Palenzuela et al. (2010a), which amounts to 
adopting Equation ( |26) with Jb = for the current and 
to Equation (|27) for ensuring the FF constraint d25). 



• driveri: denotes the first step of our "driver" approach 
and which amounts to adopting Equation ( |26) with the 
parallel component of the current specified by Equa- 
tion d29l. 



discrete 2 : denotes the second ste p of the "discrete" 
approach of Pal enzuela et al.| ( [20 1 0a) , which amounts to 
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the modification of the electric field according to Equa- 
tion ( |32) , 

• driver 2 : denotes the second step of our "driver" ap- 
proach and which amounts to adopting Equation ( (33) 
for the current. 

• continuum: denotes the continuum approach in which 
the current is specified by Equation ( (31) . 

As it is evident from Figure [T| all of the methods satisfy 
the orthogonality condition ( |2"5) essentially to machine preci- 
sion (left column). Not surprisingly, the discrete prescriptions 
discretei (combined with discrete 2 ) is particularly effi- 
cient in removing any component of the electric field parallel 
to B l , either in isolated BHs (top row) or in the case of an 
inspiralling binary (bottom row). In this latter case, the bump 
of EiB 1 at t ~ 400Af simply corresponds to the time of the 
merger and the constraint decreases after that. Similarly, the 
right column of Figure [T] shows that all prescriptions are also 
able to enforce to comparable precision the current-sheet con- 
dition of Equation ( (30) , but also that the discrete recipe ( [32) 
is slightly less effective in the case of an inspiralling binary 
(bottom right panel). 

The main conclusion to draw from Figure [T] is that, at 
least globally, all methods provide a comparable and actu- 
ally very good enforcement of the FF conditions. Their local 
performance, however, is rather different and this is shown 
in Figure [2] which reports the electrical currents as com- 
puted for a representative configuration of a single spinning 
BH with dimensionless spin a = J/Al 2 = 0.7. In the 
top panels we have reported the current vectors in the plane 
(x, y,z — 1.92 M), while in the bottom ones the currents in 
the plane (ar, y = 0, z). The two columns, on the other hand, 
contrast the currents when computed using the discretei 
and discrete 2 approaches (left column) or when computed 
using our driveri and driver 2 approaches (right column). 

A rapid comparison is sufficient to highlight that although 
both approaches yield an FF condition, the solution is very 
different, particularly on small scales. More specifically, 
when the combination of methods driver!-driver 2 is 
adopted (right column), strong meridional currents are clearly 
visible and form a jet-like structure, with negative currents in 
the central parts of the jet and positive ones on the edges of the 
jet. This current distribution is what is expected and it resem- 
bles the typical structure of the FF magnetosphere of a rotat- 
ing BH obtained through the solution o f the Grad -Shafranov 
equation (see, for instance, Figure 7 in |Beskin] ( |1997[ )). On the 
other hand, the corresponding currents when the prescriptions 
discretei and discrete 2 are used (right column) do not 
show evident signs of descending currents and, rather, they 
show unphysical features around the BH and discontinuities 
along the ~ ± 45° diagonals when seen in the (x, z) plane. 
In addition, the currents tend to be predominantly contained 
in planes which are parallel to the (x, y) plane (see the top 
row) and thus do not show the circulations which are instead 
captured with our drivers approach. 

Overall, the comparison presented in Figure |2]confirms our 
suspicions that, while providing a solution which is globally 
FF, the prescriptions discretei and discrete 2 are not guar- 
anteed to yield solutions that are locally accurate and can ac- 
tually lead to solutions with large discontinuities. For these 
reasons we believe that our approaches driveri-driver 2 
should be preferred in treatments of FF electrodynamics. As a 



final remark we also note that our prescriptions ( |29) and ( (33) 
also provide a (small) saving in computational costs. Since 
we use an algebraic prescription for the current which auto- 
matically drives the solution to the FF regime, we do not need 
to perform the expensive chec ks at every gridpo i nt that come 
with the approach suggested in Palenzuela et al. (2010a). 

6. FORCE-FREE ELECTRODYNAMICS OF BBH MERGERS 

After having discussed the details of our implementation of 
the FF conditions and having shown its higher accuracy with 
respect to alternative suggestions in the literature, in what fol- 
lows we concentrate our discussion on the FF electrodynam- 
ics accompanying the inspiral and merger of BH binaries. In 
particular, we will discuss the subtleties which emerge with 
the subtraction of the background radiation, the spatial distri- 
bution of the charge density, the EM and GW zones, and the 
scaling of the EM luminosity with frequency. 

6.1. Subtraction of Background Radiation 

As anticipated in Section [3] our measure of EM radiation 
is influenced by the choice or a uniform initial magnetic field 
within the computational domain, which leads to nonzero ini- 
tial values for $ 2 and $ - Hence, a proper identification of 
this background radiation is essential for the correct measure 
of the emitted luminosity and to characterize its properties. 

The generic expression ( |37) for the EM luminosity can 
be evaluated in combination with Equation ( f39) , that is, by 
setting as background values those of the Newman-Penrose 
scalars $2 and $0 at the initial time. Note that initial values 
of these scalars are the same they have in an ele ctro vacuum 
scenar io (they are indeed the same considered in Mosta et al. 
(2010)), and thus the "background subtraction" corresponds 
in this case to the subtraction of the EM emission coming 
from a magnetic field which is asymptotically uniform. Of 
course, the initial time is as good as any other time and we 
could in principle choose $ 2 ,b and $ ,b at any time t > 0. In 
this case, however, we would have to deal with the additional 
complication that for any choice other than t = 0, the back- 
ground radiation will also have an azimuthal modulation as a 
result of the orbital motion and hence it will not be simply an 
m = background. 

The angular distribution of the emitted radiation when pro- 
jected onto a 2-sphere, in fact, shows the presence of two jets 
but also of two extended lobes, which rotate at the same fre- 
quency of the binary and that provide the bulk of the EM 
emission (see Figure 1 of Paper I). As a result, any back- 
ground subtraction at t 7^ will also have an to = 2 com- 
ponent which will interfere with the to = 2 evolution of the 
emitted flux, introducing a modulation on the emission. The 
latter, however, will average over one orbit, leading to a net 
emitted luminosity which is the same obtained when using 
$2.b — $2(i = 0) and $ ,b = $o(t = 0). We have ver- 
ified that this is indeed the case by using background values 
at different times and obtained values of the luminosity which 
can be instantaneously different, but that once integrated over 
time yield the same emitted EM energy. As a result, the back- 
ground choice ( (39) represents by far the most convenient one. 

We have also mentioned in Section|3]that an alternative and 
equivalent estimate of the emitted EM luminosity can be ob- 
tained after removing the non-radiative parts of the emission 
(cf. expression (|40)). In order to isolate the radiative contri- 
butions from the non-radiative ones, we have reported in Fig- 
ure [3] the evolution of the real (thick lines) and of the imagi- 
nary (thin lines) parts of the I = 2, to = and I = 2, to = 2 
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FIG. 2. — Comparison of the electric currents for a single spinning BH with dimensionless spin parameter a = J/M 2 = 0.7 on the plane (x, y, z = 1.92 M) 
(top row) and on the plane (x, y = 0, z) (bottom row). All panels refer to the same time t = 102 M, when the solution has reached a stationary state. 
The currents are computed either through the fully discrete approach of dlscretei-discrete2 (left column) or through our continuous driveri-driver2 
approach (right column). While both solutions satisfy the FF condition, it is clear that the use of the drivers provides also an accurate solution. 



modes of $2 and $o- These modes are obtained from the pro- 
jection of the Faraday tensor onto the tetrad ( [35] l. Note that 
the only modes that have a regular time modulation, and are 
therefore radiative, are ($2)22 and ($0)22, while the real parts 
of the ($2)20 an d ($0)20 are essentially constant in time, indi- 
cating that these are not radiative modes, and could represent 
a way to measure the background radiation. The imaginary 
parts of ($2)20 and ($0)201 on the other hand, do show a reg- 
ular evolution in time and a ringdown, but their values are 
much smaller (i.e. , two orders of magnitude or more) and do 
not play a significant role in estimating the total radiation. 
As a result, we can write expression (|40|i explicitly as 



2,B 



Re($ 2 



$03 ^ Re($ )i 



(43) 



and in doing so we obtain an estimate that is very similar re- 
sults to the ones reported in |Neilsen et al. ( 2011| >, where ex- 
pression ( |3~8] > was used. 

As discussed in Section [3] the use of Equation p8| ) as an 
estimate of the emitted luminosity is subject to the validity of 
the assumption $2.b ~ ^o.b, or after using Equation ( |43j ), of 
Re($2)20 ~ R-e(<l>o)20- This condition is true only as a first 
rough approximation, as shown in the right panel of Figure [T 
which reports the evolution of Re(<l>2)20 (red dotted line) am 



of Re($o)20 (black dashed line), as extracted at 100 M for 
the non-spinning binary sq. Clearly, these two multipoles are 
almost constant in time and comparable, but not identical and 
their difference then affects the validity of expression ( |3~8j) . 
This consideration, together with the fact that expression (|38]> 
represents an approximation which needs to be validated a 
posteriori, leads us to the conclusion that Equation ( (37] i rep- 
resents a more accurate and robust measure of the emitted lu- 
minosity in the scenario and model considered here. 

6.2. Properties of the EM Luminosity 

Having clarified our strategy in the subtraction of the back- 
ground radiation, we present in Figure [4] a comparison of the 
evolution, measured in hours before the merger, of the lumi- 
nosities as computed with expression ([37) and either the pre- 
scriptions |39) or (j40j» for the background subtraction. 

More specifically, the thick lines refer to the total luminos- 
ity, while the thin ones to the luminosity in a polar ca p of 
5° semi-opening angle, measured using either expression ( |39] l 
(red solid line), expression ( ftO) (blue dotted line), or through 
the expression in terms of the Poynting vector ( |4l"| (black 
dashed line). The left panel refers to the binary of non- 
spinning BHs (i.e. , so), while the right one to the binary 



10 



0.001 
0.0001 

io- 5 
io- 6 
io- 7 
io- 8 
£ io- 9 

° 10-10 

o 

- io- 11 
io- ia 
io- 13 
io- 14 
io- 15 
io- 16 



Re(* 2 ) 20 
Re(* 2 ) 22 
Re(*„) 2 „ 

Re(* ) 22 



Im(* 2 ) 20 

Im(* 2 ) 22 

Im(* ) 20 

Im(* ) 22 




10" 



200 



400 
t [Ml 



600 



800 



10" 



io- 





R e(* 2 ) 20 : 




Re(* ) 20 . 


: i \ 

- : 




J i 

■) i • 

1 

1 

1 

i 1 




i 





200 



400 
t [Ml 



600 



800 
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FIG. 4. — Time evolution measured in hours before the merger of the EM luminosity at 100 M when M = 10 8 Mq and Bo = 10 4 G. The thick lines refer 
to the total luminosity, while the thin ones to the luminosity in a polar c ap o f 5° semi-opening angle, measured using either expression J35J (red solid line), 
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so), while tne right one to the binary with spinning BHs (i.e. , sg). Note that in this latter case a certain eccentricity is detectable in the EM luminosity, although 
it is much smaller in the GW luminosity. 



with spinning BHs (i.e. , s 6 ). In both cases the extraction is 
made at a distance of 100 M and the values in cgs units refer 
to a binary with a total M = 10 8 M and a magnetic field 
Bq = 10 4 G. Such magnetic-field strengths match the values 
as estimated from ra dio observations of parsec-sc ale jets in 
active galactic nuclei (O'Sullivan & Gabuzda 200"9]l. 

As expected the three measures match very well and, in par- 
ticular, the measure made with expression ([39]) is remarkably 
close to the one obtained in terms of the Poynting vector ( |4Tj i, 
that we consider the most robust measure since it involves di- 
rectly our primary evolution variables E % and B l . After the 
merger, both luminosities converge to a constant value which 
is larger than one coming from the polar-cap region (cf. thin 
lines). This is due to the fact that the background subtraction 
refers to a pure electrovacuum-condition (i.e. , uniform mag- 
netic field in a flat spacetime) and thus it does not provide an 



accurate description of an isolated spinning BH. Subtracting 
as background that of a single BH in electrovacuum would 
bring the two curves down to the values of the polar cap, but 
we have not shown this in Figure [4] to avoid a cluttering of 
curves. Note also that the measure made with expression ( |40") > 
is effectively subtracting the initial background emission and, 
at the same time, also including some incomin g radiation 
(this is true al so for the measures presented by |PalenzueIa| 



et al.| ( 2010b|a I). As a result, this measure is always (slightly) 
smaller than the one obtained with either prescriptions ( [39] l or 
( pT) . For the same reason, the contributions coming from the 
dual jets will appear comparatively larger when using Equa- 
tion ( ftO") . 

Figure H] also shows that the differences in the luminosities 
coming from the polar-cap region are instead much smaller 
and hardly noticeable. The reason behind this very good 
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agreement is simple: being integrated over a small solid an- 
gle these luminosities are not influenced by the dissimilarities 
that the different prescriptions show instead in the emitted lu- 
minosity. Overall, Figure H] shows that, as the merger takes 
place, both the diffused andthe collimated EM luminosity in- 
crease steeply, reaching values at the merger which are about 
50 times larger than the corresponding ones a few orbits be- 
fore the merger. The growth in the diffused luminosity, how- 
ever, is larger than the one in the collimated luminosity and 
the difference in the two, which was already present at the be- 
ginning of the simulations, increases as the inspiral proceeds. 
As a result, at the merger the non-collimated (total) emission 
is ~ 100 times larger than the collimated one, reaching values 
L EM ~ 10 45 ergs" 1 for a 10 8 M binary^ 

A few comments should be reserved about the different spa- 
tial distributions of the EM fluxes that come with the different 
prescriptions for the subtraction of the background radiation 
and that are erased when computing the luminosities as inte- 
gral quantities. First of all we note that the EM flux in Equa- 
tion ([37) is not necessarily positive on the 2-sphere and that 
(small) negative contributions can appear (see Figure 1 of Pa- 
per I and the corresponding color bar). These emissions, how- 
ever, do not represent a radiative field and average to zero over 
one or bit (this point was already remarked in Palenzuela et al. 
(2010), where a toy model within the membrane paradigm 
was used for the binary). This non-radiative part is far from 
being uninteresting as it could lead to a different secondary 
emission as the EM fields interact with the plasma. Unfortu- 
nately, by construction, it is impossible to investigate such an 
emission within our FF approach, but this is clearly an aspect 
of this research that deserves further investigation. Second, 
as already remarked in Paper I, while the EM fluxes do con- 
tain a dual-jet structure and even if the fluxes at the jets are 
^8 — 2 times larger than elsewhere, the global spatial dis- 
tribution is effectively dominated by a non-collimated emis- 
sion of quadrupolar nature, drastically changing the prospects 
of the detectability of the dual jets (see also discussion be- 
low). Finally, the local EM flux from the jets can in princi- 
ple be enhanced if the BHs are spinning and, indeed, within 
a Blandford-Znajek process one expects that the luminos- 
ity fro m the jets increases quadratically with the spin of the 
BH ( |Blandford & Znajek|1977||Palenzuela et al.|2010b) . The 
differences introduced by the spin are reported in the right 
panel of Figure HI which refers to the binary sq and thus with 
BHs having a dimensionless spin of J/M 2 ~ 0.6. Clearly, 
both the collimated and the non-collimated emission show a 
behavior which is similar to the one seen for the sq binary, 
with only a 50% enhancement of the EM radiation, both in the 
total and in the collimated emission (note that the two panels 
in Figure |4] have the scale). This result is the consequence 
of the fact that most of the radiation that is produced is dif- 
fused and produced by the interaction between the BH orbital 
motion and the background magnetic field. Indeed, we find 
th at the emission in th e electrovacuum evolution as computed 
Mosta et aT7| ( |2Q10| > is co mparable to the FF one (this is dif- 



m 



ferent from what reported in|Palenzuela et al." ( 2010b|a|>). The 



local spin enhancement in the dual jets is therefore present, 
but still much smaller than the diffused emission, which re- 
mains the predominant one at these separations. 
It is always useful to remark that by construction the 

8 Note that the local flux of the collimated emission can be ~ 8 — 2 times 
larger than the one in the diffused emission. However, being limited to a very 
small solid angle, the corresponding luminosity is 100 times smaller. 
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FIG. 5. — Evolution of the EM (top panel) and the GW luminosity (bottom 
panel) integrated over 2-spheres located respectively at r = 20 , 100, and 
180 M. Thick lines refer to the diffused emission, while thin ones to the 
emission from a polar cap of 5° semi-opening angle; the data refer to the 
spinning se, binary and both the EM and the GW luminosities are computed 
including modes up to the I = 8 multipole. Note that the gravitational- 
wave zone is already well defined at 100 M, while the EM one is not even at 
180 M. 



Newman-Penrose scalars, either for the gravitational sector, 
i.e. , ^4, or for the EM one, i.e. , $0,^2, provide non- 
ambiguous quantities only at very large distances from the 
sources, that is, in the corresponding "wave zone". It is obvi- 
ous then that any measure of such radiation quantities in the 
strong-field region, risks to be incorrect. Less obvious is how- 
ever the fact that the wave zones can be different whether one 
is considering the gravitational or the EM radiation, with the 
latter starting at considerably larger distances than the former. 
This is summarized in Figure B] which reports the EM (top 
panel) and the GW luminosity (bottom panel) integrated over 
2-spheres located respectively at r = 20, 100, and 180 M. 
The data refer to the spinning sq binary, with both the EM 
and the GW luminosities having been computed including 
modes up to the I — 8 multipole; thick lines refer to the dif- 
fused emission, while thin ones to the emission from a polar 
cap of 5° semi-opening angle. Clearly, the estimates made at 
r = 20 M in both channels are rather different (and incor- 
rect) from those made at larger radii, where the radiation has 
reached its wave-like solution. Also striking is that while the 
GW estimates at 100 M and 180 M are essentially indistin- 
guishable (bottom panel), the corresponding ones in the EM 
channel are not yet identical. This indicates first that the GW 
zone is much closer than the EM one and reached already 
at r ~ 100 M, and, second, that extraction radii larger than 
r ~ 200 M should be considered when measuring the EM 
radiation. We note that the evidence of a relative "proximity" 
of the GW zone to the strong-field dynamical region of space- 
time is somewhat surprising, but also in substantial agreement 
with the bulk of evidence emerging in favor of a description of 
the dynamics of the BHs which is very well described by PN 
or other approximation techniques. This good agreement is 
indeed perfectly understandable if the weak-field wave zone 
starts only a few tens of M away from the BHs. 

6.3. Frequency Scaling 
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FIG. 6. — Left Panel: frequency scaling for the non-spinning binary so of the GW luminosity rescaled of a factor 10 (black solid line), of the diffused 
EM luminosity (red solid line), and of the collimated EM luminosity computed in a polar cap with a semi-opening angle of 5° (blue solid line). Note that the 
diffused EM luminosity has a behavior which is compatible with Q 1( V3— 8/3 as (j oes as me luminosity. The collimated EM luminosity, on the other hand 
has a scaling compatible with £7 5 / 3- 6 / 3 . Right Panel: the same as in the left panel but reporting only the GW emission and extrapolating back in the past to 
determine when the collimated and the diffused emissions are comparable. For a binary with 10 s Mq this happens ~ 21 days before merger. 



As remarked already in Paper I, an accurate measure of 
the evolution of the collimated and non-collimated contribu- 
tions of the emitted energies is crucial to predict the proper- 
ties of the system when the two BHs are widely separated. 
This measure, however, is all but trivial as it requires a re- 
liable disentanglement of the collimated emission from the 
non-collimated one and from the background. We have seen 
in Figure|4]how the total EM luminosities show a very similar 
evolution as long as sensible subtractions of the background 
radiation are used. We have also discussed that independently 
of the choice made, the diffused emission is mostly quadrupo- 
lar and hence with a dependence that is the same as the GW 
ft 10 / 3 



Palenzuela et al. 



one, i .e. , ~ il^'^, as a l ready shown by 
( |2010[ ) and Mosta et al] ( |2010j l. Figure |6j considers more 
closely this issue by reporting in the left panel the change 
of the different gravitational and EM luminosities in the or- 
bital evolution as a function of the GW frequency ft GVJ . More 
specifically, we report the diffused EM radiation as computed 
with expressions ( [37] i and ( |40| (red solid line) and the colli- 
mated emission when computed over a polar cap with a semi- 
opening angle of 5° (blue solid line). Also shown is the evo- 
lution of the GW luminosity (black solid line) scaled down of 
a factor 10~ 10 to make it comparable with the other luminosi- 
ties (we recall that the efficiency in GW emission is ~ 13 o r- 
ders of magnitude larger as first shown in Mosta et aL] ( |2010] l). 
The short-dashed, dotted, and long-dashed lines show instead 
the different scalings (note the figure is a log-log plot). 

It is then straightforward to realize that at the separations 
considered here the diffused emission shows a scaling with 
frequency which is L^° n ~ co11 w ^10/3-8/3^ t ^ us compatible 
with the scaling shown by the GW emission. The collimated 
emission, however, has a slower growth, with a scaling that is 
L™J' w fi 5 / 3 ~ 6 / 3 . This is different from the predicted scaling 
of L coU s» fi 2 / 3 suggested in 



(|2010b|, and 
inis differ- 



Palenzuela et al 
that we show with a light-blue long-aasned line.' 
ence is probably due to the fact that the estimate in Palenzuela 



|et al.| ( f2"010b| l was made by studying the behavior of boosted 
BHs and then extrapolating the result to the case of orbiting 
BHs. The scaling ~ f2 2 / 3 is clearly incompatible with our 
data and we suspect the accelerated motion of the BHs to be 
behind this difference and longer simulations will be useful to 
draw robust conclusions. 

Given that the diffused and the collimated emissions scale 
differently with frequency and using the rough estimates 
made above for their scaling at earlier times, [j we can de- 
termine the frequency (or time) when the collimated emission 
will be dominant relative to the diffused one. This is shown 
in the right panel of Figure |6| which is the same as the left 
one but where we extrapolate the scaling back in frequency. 
Our rough estimate is therefore that the collimated emission 
will be larger than the diffused one at an orbital frequency 
f2 = ^J1 GW ~ 3.2 x 1CT 5 Hz and thus ~ 21 days before the 
merger. If the conditions are optimal and the binary is oriented 
in such a way that the dual-jet system points toward the Earth, 
the luminosity from the binary would therefore be modulated 
on timescales t < 1/51 ~ 8.6 hr and smaller. While this is an 
exciting possibility, we should also bear in mind that, when 
extrapolated back to the time when it becomes dominant, the 
collimated emission has also decreased by almost one order 
of magnitude and to luminosities that are only of the order of 
~ 10 ergs -1 . Luminosities ~ 10 45 ergs _1 are also typ- 
ical of radio-loud galaxies and thus the determination of an 
EM counterpart can be challenging if such sources are near 
the candidate event. Clearly, the bottom line of these con- 
siderations is that longer simulations need to be performed to 
assess the early-inspiral scaling of the different luminosities 
and more realistic scenarios need to be considered to assess 
whether the collimated or the diffused emission can serve as 



9 In reality we expect the scaling with frequency to be different in the 
different stages of the inspiral, just as it is the case for the GW emission. 
However, as a first approximation we can assume that the frequency does not 
change significantly in the early stages of the inspiral. 
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FIG. 7. — Small-scale two-dimensional distribution of the charge density for a sq binary in the early inspiral phase at t = 89 M (left column), at merger 
t = 672 M (middle column), and at ringdown t = 800 M (right column). The top panels show the charge density in the (x, y) plane, while the bottom ones in 
the (x, z) plane. Visualizations artifacts appear as thin stripes at the boundaries between refinement levels; the data in those stripes are of course regular. 



an EM counterpart t o the merger o f binary syste m of super- 
massiv e BHs (see jGiacomazz o et al.| ( |20T2] l and [Noble et al.| 
(2012) for some recent attempts). 

6.4. Charge-density Distribution 

In this concluding section we concentrate on the spatial dis- 
tribution of the charge density produced during the inspiral 
and merger, providing informa tion which is complementary to 
the one already presented by Palenzuel a et al. (2010); Palen- 
|zuela et al.| ( |2010b| > and |Neilsen et al.| ( |201 1) . We recall that 
in our simulations the charge density is not an evolutionary 
quantity, but, rather, it is computed from the constraint equa- 
tion (l7j. We also recall that because we are very effective in 
enforcing the FF condition (see discussion in Section B), we 
cannot fully explore the physical consequences of the charge 
distribution we produce. This is because in the most inter- 
esting regions of these distributions, that is, in those regions 
with no (or very small) net charges and which are remi niscent 
of the vacuum-gap regions in pulsar magnetospheres (Becker 
2009), the electric field along the magnetic field will be zero 
to machine precision and hence it will not be able to accelerate 
particles to very high Lorentz factors (as instead is expected 
in the polar regions of pulsar magnetospheres). To further 
limit the amount of information that can be extracted directly 
from our simulation is the fact that an FF code does not allow 
for an unambiguous calculation of the plasma velocity, which 
can only be estimated a posteriori based on a certain n umber 
of assumptions. As an example, Hirotani & Oka moto1 ( |1998| l 
argued that it is possible to compute the final Lorentz factor of 
a plasma in an FF magnetosphere if there is a non-negligible 
component of the parallel electric field and a radiation drag 
dominated by Thompson scattering. 

In spite of these limitations, the charge-density distribution 
remains a very interesting quantity and we have reported it in 



Figures [7jand[8] The three top panels of Figure[7J in particular, 
show the charge distribution on the (x, y) plane, while the bot- 
tom ones on the (a;, z) planes at three different instants in the 
evolution of the spinning binary ,s 6 . More specifically, in the 
early inspiral phase (t = 89 M), at the merger (t = 672 M), 
and at ringdown (t = 800 M). The color code highlights the 
presence of positive (red) and negative (blue) charges, which 
are produced both because of the orbital motion of the BHs, 
but also because of the intrinsic spin of the BH. The first con- 
tribution can be appreciated from the first two columns of Fig- 
ure 17] while the second contribution is the only one responsi- 
ble for the charge distribution in the last column. Much of 
this distribution of charg es can be easily int erpreted within 
the membrane paradigm ( |Thorne et al.| 1986] > as the result of 
an effective Hall effect arising when the BH horizon (i.e. , 
the "membrane") moves, either as a result the orbital motion 
or through its spinning motion, across a magnetic field. In 
analogy with the classical Hall effect, a charge separation will 
be produced as shown i n Figure [7] (see a lso the discussion 
in |Neilsen et"aT1 ( |2011| i; |Lyutikov| ( |2011| i). Note that since 
they both refer to isolated spinning BHs (although with differ- 
ent spins), the right column of Figure [7J should be compared 
with the right column of Figure [2] which shows instead the 
electric currents. 

Additional information is shown in Figure IS] where the 
charge-density distribution is rendered in three dimensions at 
the same representative times shown in the panels of Figure|7j 
and on much larger length scales. This representation high- 
lights that the distribution is far more complex than a sim- 
ple dual-jet structure and is instead typical of a double-helical 
symmetry, simil ar to the p attern for the Poynting flux shown 
in Palenzuela et al. (2010b a). Although it is not possible to 
investigate further, within an FF approach, the consequences 
of this regular and alternate distribution of positive and neg- 
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FIG. 8. — Top row: large-scale three-dimensional distribution of the charge density for the se, binary in the early inspiral phase at t = 89 M (left panel), at the 
merger t = 672 M (middle panel) and at ringdown t = 800 M (right panel). In these panels only the largest values of the charge density are shown. Bottom 
row: three-dimensional distribution of the charge density at ringdown only, t = 800 M. Starting from the left, the panels show smaller and smaller values 
of the charge density, revealing a much more extended conical-shaped structure with a double-helical distribution of opposite charges. Clearly, charge-density 
distribution is far more complex than what would be deduced from the top panels only. 



ative charges, it is clear that it can lead to rather intriguing 
particle acceleration processes along the surfaces separating 
regions of different charges. The resulting accelerated parti- 
cles could further cascade into less energetic charges and lead 
to a potentially detectable emission. 

It is worth remarking, however, that the charge-density dis- 
tribution is not restricted to a small cylindrical area compris- 
ing the two inspiralling BHs, as it may erroneously appear 
from the top panels of Figure [8] and which shows only the 
regions where the charge density is the largest. Rather, it in- 
volves the whole region in causal contact with the binary, as 
shown in the lower panels of Figure [8] which refer instead to 
the ringdown phase only (t = 800 M). Starting from the left, 
the different panels are drawn exhibiting increasingly smaller 
values of the charge density and thus revealing a much more 
extended conical-shaped structure with a double-helical dis- 
tribution of opposite charges at its core. Additional investiga- 
tions away from the FF regime will be necessary to assess the 
astrophysical impact of these structures. 

7. PROSPECTS AND CONCLUSIONS 

Assessing the detectability of the EM emission from merg- 
ing BH binaries is much more than an academic exercise. The 
detection of EM counterpart, in fact, will not only act as a 
confirmation of the GW detection, but it will also provide a 
new t ool for testing a numb er of fundamental astrophysical 
issues Haiman et al. ( 2009| l. In particular, it will offer the 
possibility of testing models of galaxy mergers and accretion 
disks, of probing basic aspects of gravitational physics, and 



of dete rmining c osmological parameters once the redshift is 
known (Phinney 200~9jl. 

Computing reliable estimates from this scenario is made 
difficult by the scarce knowledge of the physical conditions 
in the vicinity of the binary when this is about to merge. Nev- 
ertheless, relying on a number of assumptions with varying 
degree of realism, several investigations have been recently 
carried out to investigate the properties of these EM coun- 
terparts either during the stages that precede the merger or 
in those following it. As an example, several authors have 
recently conside red the interaction between t he binary and a 
dense gas cloud ( Armitage & Natarajan 2002; van Met er et al.| 
20T0l|Bode et al'|2010||Farris et al.|2010||Lodato et al.|2009| 
Chang et al.||2010l IFarris et al.||2011MBode et al.||2012| |Gi-| 
acoma zzo et al.|20l2| [Noble et al.||20l2p even though astro- 
physical considerations seem to suggest that during the very 
final stages of the merger the SMBBH will inspiral in a rather 
tenuous intergalactic medium. At the same time, scenarios 
which do not involve dense matter distributions in the vicin- 
ity of the binary have also been considered. In these cases, the 
SMBBH is assumed to be inspiralling in electrovacuum and in 
the presence of an extern al magnetic field which is anchored 
to the circumbinary disk (Palenzuel a et al.|20 09; Most a et aL] 
2010| > and the energy emitted in EM waves is ~ 13 orders of 



magnitude smaller than the one emitted in GW for a typical 
binary of supermassive BHs wit h mass M = 10 8 M in an 
ambient magnetic field of 10 4 G ( |Mosta et al.|2010| >. 

Furthermore, when charges and current s are considered 
within an FF regime, the numerical results of Pale nzuela et al!] 
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(2010b a) have shown that, if taking place in a uniform mag- 
netic field, the merger event would be accompanied by the 
EM emission from a dual-jet structure, acting as a finger- 
pri nt of the me rger itself. A detailed analysis carried out 
in Kaplan et al. ( 2011| l addressed the problem of whether 
such merger flares can be detected by ongoing and planned 
wide-field radio surveys, su ch as the Square Kilometer Array 
pathfinder (Johnsto n et al.|[2007 i. The conclusion was that, 
owing to the short timescales associated with the merger, no 
more than one event per year would be det ectable by such 
blind surveys. In a recent paper (Moesta et al.|2012]l we have 
revisited the estimates made in Palenzuela et al.| ( |2010b|a| l 
and shown that while a dual-jet structure is present during 
the inspiral, and while the fluxes can be larger near the jet, 
the collimated luminosity is subdominant of a factor ~ 100 
with respect to the total luminosity, which is instead predom- 
inantly quadrupolar. Furthermore, spin-related enhancements 
are only very small and less than 50% when considering a 
spinning binary with dimensionless spins J/M 2 = 0.6. 

Our results have been obtained adopting a consistent mea- 
surement of the EM luminosity and an improved numerical 
strategy for the treatment of the FF condition, both of which 
have been discussed in detail in this paper. More specifically, 
we have shown that we do not implement the FF condition at 
a discrete level, but rather we obtain it via a damping scheme 
which drives the solution to satisfy the correct condition. This 
difference is important for a correct and accurate description 
of the current sheets that can develop in the course of the 
simulation. We have also studied in greater detail the three- 
dimensional charge distribution produced as a consequence 
of the inspiral and shown that it possesses a complex but or- 
dered structure with a double-helical distribution of opposite 
charges tracing the motion of the two BHs. 

Although our simulations show that the dual-jet structure 
is subdominant on the timescale over which the simulations 
have been carried out, they also indicate that the growth rates 
of the collimated and diffused luminosities are different, thus 
suggesting that sufficiently early in the inspiral the collimated 
emission will be the dominant one. Computing accurately 
these scaling rates is of course crucial since it allows for the 
determination of the time during the inspiral in which the dual 
jets are dominant could modulate the emission if the binary is 
suitably oriented. When consi dering the observa tional impli- 
cations of this possibility, O'Shaughnessy et al. ( 201 l|l have 
concluded t hat future blind radio surveys like VAST (Banyer 



|et al.|[20T2] > would easily detect the effects of these modula- 
tions, with a frequency of up to one per day. 

We have therefore provided the first quantitative estimates 
of the scaling of the EM emission with frequency and shown 
that the diffused part has a dependence that is very close to 
the one exhibited by the GW luminosity and therefore of the 
type L non - co11 « The collimated EM emission, 

J 1 EM ' 



on the other hand, scales like L coU « SI 5 / 3 6 / 3 , thus with 

' EM ' 

a steeper dependence than L™' 1 » il 2 / 3 , as previously sug- 
gested by Palenzuela et al.| ( |2010b |. In light of these scalings 
and considering a non- spinning binary, we conclude that the 
collimated emission will be larger than the diffused one at an 
orbital frequency of ~ 3.2 x 10~ 5 Hz and thus ~ 21 days 
before the mergerp*] When this happens, the collimated lumi- 
nosity will be about an order of magnitude smaller than the 
one considered here and of the order of ~ 10 42 ergs _1 for 
a typical 10 8 M & binary in a magnetic field of 10 4 G. Such 
a luminosity is about 1000 times smaller than the typical lu- 
minosity of radio-loud galaxies and thus determination of an 
EM counterpart can be challenging if such sources are near 
the candidate event. 

As a concluding remark we note that while our study ad- 
dresses several points which were not fully investigated be- 
fore, it also leaves open a number of questions. One of 
these questions is the efficiency of the secondary emission 
that could be generated either by the diffused component or 
by the collimated one. The richly complex structure of the 
charge-density distribution, in fact, can be the site where even 
small electric fields along the magnetic field lines would be 
able to accelerate particles to very high Lorentz factors, lead- 
ing to a secondary emission similar to the one expected in the 
polar regions of pulsar magnetospheres. Unfortunately, how- 
ever, our use of an FF condition (and our ability to maintain it 
essentially to machine precision) prevents us from producing 
such electric fields and hence the corresponding accelerations. 
Another and related unresolved issue is the fate of the Poynt- 
ing flux once it impacts the intergalactic medium. Even in the 
optimistic case in which the majority of the Poynting flux is 
converted into radio emission via synchrotron processes, the 
EM radiation (either collimated or diffused) will eventually 
exit the evacuated central region around the binary and pene- 
trate in the ambient medium. When this happens, part of the 
Poynting flux will be converted into kinetic energy and repro- 
cessed in several EM wavebands, not necessarily in the radio 
range. Clearly, Ion ger simulations and more realistic sce- 
narios are needed to shed further light on the properties of the 
EM counterpart to the inspiral and merger of binary of super- 
massive BHs. 
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APPENDIX 

ON THE IMPLEMENTATION OF THE IMEX SCHEME 

The prototype of the stiff system of partial differential equations can be written as 

d t U = F (U) + aR(U) , 



(Al) 



10 Clearly, this equivalence in the emission will take place much earlier 
(and at smaller luminosities) if the scaling is less steep than ~ f2 10 / 3 . 



nuclei suggest that in these cases more than 70% of the Poynting flux can 
be converted into kinetic e nergy leading to flows with Lorentz factors of the 



1 ' Numerical MHD simulations in the context of jets from active galactic order of V ~ 10 jKomissarov et al.|2007| 
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where 1/a > is the relaxation time. In the limit (T4oo the system becomes stiff, since the relaxation of the stiff term R(U) 
is very different from the timescale of the non-stiff part F(U). 

The evolution of the electric field ( fl7| i becomes stiff for high values of the conductivity as in the Ohm law ( p6*| ). We perform 
a split of its right-hand side in potentially stiff terms and regular ones, 

d t E = F E + R E , (A2) 

where 

F E = e^ k e 4 * [ (dj a ) % k B c + a ( 4 % k 8 3 <j> + dj % k ) B c 

^ k E,B k 



a 7 c fe dj B c }+ CpE 1 - aK E l - aq- 



B 2 



R E = -aJ B —. (A3) 

A solution for the magnetic field is obtained by evolving Equation ( p"8j ) using only the explicit part of the Runge-Kutta solver. 
The evolution of the electric field uses both the explicit part of the Runge-Kutta solver (see Table 1) for the F E and the implicit 
part for R E (see Table 2), and leads to an approximate solution {-E*}. The full solution requires inverting the implicit equation 

E = E* + a u AtR E (E), (A4) 

which depends on the fields {B, E*}. 
In the case of the Ohm law ( |2"9| l the stiff part is linear in E, so an analytic inversion can be performed 

E i =(M k i )- 1 Ei, (A5) 

B i 



M k l = 5 k l + an At a a B B k . (A6) 

However, in the case of the Ohm law f33| l, the inversion is more involved as the stiff part is not linear in E. We use the 
following simplified inversion: 

E i = (M k i )- 1 E*, (A7) 

M k l = 4 4 + an At a a B \ B k |^ + 5 k \E% - B 2 ) || 

In the above equations, At is the timestep and are the diagonal coefficients of the implicit part of the RKIMEX matrix, 
whose tableau for the explicit and explicit-implicit IMEX-SSP3(4,3,3) L-stable scheme are reported below: 
where 

a = 0.24169426078821 , (3 = 0.06042356519705 , 
77 = 0.12915286960590. 



TABLE 1 

Explicit IMEX-SSP3(4,3,3) L-stable Scheme 
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TABLE 2 

Implicit IMEX-SSP3(4,3,3) L-stable Scheme 
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